###########################csc###################
#   Replication for 
#
# 	The Risks of Election Observation: International Condemnation and Post-Election Violence
# 	Inken von Borzyskowski
#   International Studies Quarterley
# 
#	  Created: 11 December 2018
#
#   Code for Figure 2, Figure A1, Table A15, Figure A3, Figure A4
#
##############################################

# install if you have not yet:
# install.packages("Matching", dependencies=TRUE) 
# install.packages("readstata13")

#Set working directory: setwd()

#######################
# Figure 2
#######################
rm(list=ls(all=TRUE))
library(foreign)

# Simulation estimates (see Stata code):                  coef             SE      CI lower    CI upper
# NegReport 0 SeriousFraud 0 -> Pr(PostElViolence=1) |   .0158698     .0188866     .0015626    .0683706
# NegReport 1 SeriousFraud 0 -> Pr(PostElViolence=1) |   .1909332     .1238674     .0352637    .4998722
# NegReport 0 SeriousFraud 1 -> Pr(PostElViolence=1) |   .0681363     .0444732     .0166228    .1908475
# NegReport 1 SeriousFraud 1 -> Pr(PostElViolence=1) |   .5348327      .158626     .2045016    .8123026

xrange <- c(.5,7.5)                                        
yrange <- c(0,1)     

png(file="Figure2.png", res=300, height=1700, width=2700)  
par(mar=c(4,4,2,2)+.1)
plot(xrange, yrange, type="n", xlab="",
     ylab="Probability of Post-Election Violence",yaxp=c(0,1,5), las=1, xaxt="none", ps=6,  main="") 
axis(1, at=1:7, labels=c("No Condemnation \n No serious fraud", "","Condemnation \n No serious fraud", "","No Condemnation \n Serious fraud", "", "Condemnation \n Serious fraud"), line=2, tick=FALSE, cex.axis =1)    
abline(h =0.0, untf = FALSE,lty="dashed", lwd=.9, col="darkgray")
abline(h =0.2, untf = FALSE,lty="dashed", lwd=.9, col="darkgray")
abline(h =0.4, untf = FALSE,lty="dashed", lwd=.9, col="darkgray")
abline(h =0.6, untf = FALSE,lty="dashed", lwd=.9, col="darkgray")
abline(h =0.8, untf = FALSE,lty="dashed", lwd=.9, col="darkgray")
abline(h =1,   untf = FALSE,lty="dashed", lwd=.9, col="darkgray")
lines(c(1,3), c(.0158698, .1909332), type="l", lwd=2, col="darkgray", pch=20, cex=.5) 
lines(c(5,7), c(.0681363, .5348327), type="l", lwd=2, col="darkgray", pch=20, cex=.5) 
text(2,0.16,label= "+ 0.17")
text(5.9,0.37,label= "+ 0.46")
lines(c(1,1),c(.0015626, .0683706), type="l", cex=1,  lwd=3, col="black")   
lines(c(3,3),c(.0352637, .4998722), type="l", cex=1,  lwd=3, col="black")     
lines(c(5,5),c(.0166228, .1908475), type="l", cex=1,  lwd=3, col="black")  
lines(c(7,7),c(.2045016, .8123026), type="l", cex=1,  lwd=3, col="black")    
points(c(1,3,5,7),c(.0158698, .1909332, .0681363, .5348327), pch=20, cex=1, lwd=3, col="black")              
points(c(1,3,5,7),c(.0015626, .0352637, .0166228, .2045016), pch="--", cex=2, lwd=5, col="black")                                                   
points(c(1,3,5,7),c(.0683706, .4998722, .1908475, .8123026), pch="--", cex=2, lwd=5, col="black")  
text(0.7,0.03,label= "0.02")
text(3.3,0.18,label= "0.19")
text(4.7,0.07,label= "0.07")
text(7.3,0.53,label= "0.53")
dev.off()



#######################
# Figure A1
#######################
rm(list=ls(all=TRUE))
library(foreign)
dat<-read.csv("FigureA1data.csv")
attach(dat)

png(file="FigureA1.png", res=300, height=1200, width=2500)  
par(mar=c(4,4,1,1)+.1,xaxs="i", yaxs="i")
hist(dat$GapElReport, xlim=c(-6,35),ylim=c(0,30), col="red", breaks=seq(-6,35,l=42), xlab="Number of Days after Election (0 = Election Day)",ylab="Number of Elections",main="", xaxt="n",yaxt="n")
hist(dat$GapElResult, add=T, col=rgb(0, 0, 1, 0.3),breaks=seq(-6,35,l=42),  xaxt="n",yaxt="n")
axis(1, at=c(-5.5,-0.5,4.5,9.5,14.5,19.5,24.5,29.5,34.5,39.5), labels=c("-5","0","5","10","15","20","25","30","35","40")) 
axis(2, las=1,at=c(0,10,20,30), labels=c("0", "10","20", "30")) 
legend(20,25, c("IO Report", "Election Result"), cex=1.1,fill=c("red", rgb(0, 0, 1, 0.3)),bty="n")
box()
dev.off()
  
  

#######################
# Table A15 and Figure A3 
#######################
rm(list=ls(all=TRUE))
library(foreign)
library(xtable)
library(readstata13)
library("Matching")

dat<-read.dta13("Condemnation.dta")
attach(dat)
dat1 <- subset(dat,dat$observedW46 ==1, select = c(PostElViolence, NegReport, SeriousFraud, LoserVoteShare, Boycott, 
        type3, previolev, DPI_yrsoffc_Lag1Log, Prio_CWPast10, WB_GDPpc_Lag1Log, Ross_oilgas_valuePOP2000_Lag1Log, 
        P_xconst_Lag1, WB_PopM_Lag1Log, Finkel_EF_Lag1, WB_ODAPCconstant_Lag1Log, FirstMultiParty, Lindberg_Turnover, 
        COUNTRY_NAME, year))
dat1<-na.omit(dat1)
attach(dat1)
rm(dat)
row.names(dat1) <- NULL
dat1$row.names = 1:nrow(dat1)

# get propensity score
mod <- glm(NegReport ~ previolev + SeriousFraud + WB_GDPpc_Lag1Log +
             WB_ODAPCconstant_Lag1Log + FirstMultiParty + as.numeric(ordered(Lindberg_Turnover)), 
           family=binomial(link="logit"), data=dat1)
summary(mod)
dat1$propensity  <- mod$fitted

set.seed(54321)

# set DV, treatment, balance variables 
Y  <- PostElViolence
Tr  <- NegReport
X = cbind(dat1$propensity, SeriousFraud, LoserVoteShare, Boycott, +
            type3, previolev, DPI_yrsoffc_Lag1Log, Prio_CWPast10, WB_GDPpc_Lag1Log,  + 
            Ross_oilgas_valuePOP2000_Lag1Log, P_xconst_Lag1, WB_PopM_Lag1Log, Finkel_EF_Lag1)
BalanceMat <- cbind(dat1$propensity, SeriousFraud, LoserVoteShare, Boycott, +
            type3, previolev, DPI_yrsoffc_Lag1Log, Prio_CWPast10, WB_GDPpc_Lag1Log,  + 
            Ross_oilgas_valuePOP2000_Lag1Log, P_xconst_Lag1, WB_PopM_Lag1Log, Finkel_EF_Lag1)

# genetic matching
genout <- GenMatch(Tr=Tr, X=X, BalanceMatrix=BalanceMat, estimand="ATT", M=1,
            pop.size=2000, max.generations=10, wait.generations=1)

# get result/estimate ATT (referenced on Appendix page 25) 
rr <- Match(Y=Y, Tr=Tr, X=X, estimand="ATT", Weight.matrix=genout)
summary(rr)

# check balance
mb <- MatchBalance(NegReport~dat1$propensity + SeriousFraud + LoserVoteShare + Boycott  +
            type3 + previolev + DPI_yrsoffc_Lag1Log + Prio_CWPast10 + WB_GDPpc_Lag1Log + 
            Ross_oilgas_valuePOP2000_Lag1Log + P_xconst_Lag1  + WB_PopM_Lag1Log + Finkel_EF_Lag1,
            match.out=rr, nboots=500)

# save balance estimates
I <- 13
names <- c("propensity", "SeriousFraud", "LoserVoteShare", "Boycott",
           "type3", "previolev", "DPI_yrsoffc_Lag1Log", "Prio_CWPast10", "WB_GDPpc_Lag1Log",
           "Ross_oilgas_valuePOP2000_Lag1Log","P_xconst_Lag1","WB_PopM_Lag1Log","Finkel_EF_Lag1")
before.mat <- matrix(NA,13,5)
after.mat <- matrix(NA,13,5)
before <- mb$BeforeMatching

for(i in 1:I){
  before.mat[i,1] <- names[i]  
  before.mat[i,2] <-1
  before.mat[i,3] <- mb$BeforeMatching[[i]]$qqsummary$meandiff
  before.mat[i,4] <- ifelse(is.null(mb$BeforeMatching[[i]]$ks$ks.boot.pvalue),mb$BeforeMatching[[i]]$tt$p.value,mb$BeforeMatching[[i]]$ks$ks.boot.pvalue)
  before.mat[i,5] <- ifelse(before.mat[i,4]<=0.05,1,0)
  after.mat[i,1] <- names[i] 
  after.mat[i,2] <-2
  after.mat[i,3] <- mb$AfterMatching[[i]]$qqsummary$meandiff
  after.mat[i,4] <- ifelse(is.null(mb$AfterMatching[[i]]$ks$ks.boot.pvalue),mb$AfterMatching[[i]]$tt$p.value,mb$AfterMatching[[i]]$ks$ks.boot.pvalue)
  after.mat[i,5] <- ifelse(after.mat[i,4]<=0.05,1,0)
}

balance.out <- rbind(before.mat, after.mat)
colnames(balance.out) <- c("","stage","mean_eCDF_diff","KS Bootstrap p-value","significant")

write.csv(balance.out, file = "FigureA3data.csv")
dat<-read.csv("FigureA3data.csv")
dat$KS.Bootstrap.p.value<-format(dat$KS.Bootstrap.p.value, scientific = FALSE)
dat$significant<-ifelse(dat$KS.Bootstrap.p.value<=0.05,1,0)
write.csv(dat, file = "FigureA3data.csv")
## use these data for Figure A3 below


## Table A15
# make dataframe for treated and control
rr<-data.frame(rr$index.treated,rr$index.control)
colnames(rr)<-c("treated","controls")
dat1$rowID<-as.numeric(row.names(dat1))

# recode poll type from numeric into Legislative, Executive, General
dat1$type3[dat1$type3==0] <- "L"
dat1$type3[dat1$type3==1] <- "E"
dat1$type3[dat1$type3==2] <- "G"

# join treated and control ID to original characteristics
rr$TreatedCountry <- dat1$COUNTRY_NAME[match(rr$treated, dat1$rowID)]
rr$ControlCountry<-dat1$COUNTRY_NAME[match(rr$controls, dat1$rowID)]
rr$TreatedYear <- dat1$year[match(rr$treated, dat1$rowID)]
rr$ControlYear<-dat1$year[match(rr$controls, dat1$rowID)]
rr$TreatedPolltype <- dat1$type3[match(rr$treated, dat1$rowID)]
rr$ControlPolltype<-dat1$type3[match(rr$controls, dat1$rowID)]
rr
rr <- subset(rr, select = c(TreatedCountry, TreatedYear, TreatedPolltype, ControlCountry, ControlYear, ControlPolltype))
print(xtable(rr, type = "latex", digits=0), include.rownames=FALSE, file = "AppendixTableA15.tex")



#######################
# Figure A3
#######################
rm(list=ls(all=TRUE))
library(foreign)
dat<-read.csv("FigureA3data.csv")

pdf(file = "FigureA3.pdf",
    width = 6, height = 6, onefile = TRUE, family = "Times")
par(mar=c(4,4,2,2)+.1)
plot(dat$stage,dat$mean_eCDF_diff, ylim = c(0, 1), xlim = c(0, 2.5), 
     ylab = "Standardized Difference in Means", xaxt="n", 
     xlab="", pch=21, 
     bg=ifelse(dat$sig==1,
               "black", "gray"), 
     col=ifelse(dat$sig==1,
                "black", "gray"))
axis(side = 1, at = c(1,2), labels=c("All Data", "Matched Data"))
lines(c(1.03,1.97), c(dat$mean_eCDF_diff[dat$X==1], dat$mean_eCDF_diff[dat$X==14]), type="l", lwd=.8, col="darkgray", pch=20, cex=.5) 
lines(c(1.03,1.97), c(dat$mean_eCDF_diff[dat$X==2], dat$mean_eCDF_diff[dat$X==15]), type="l", lwd=.8, col="darkgray", pch=20, cex=.5) 
lines(c(1.03,1.97), c(dat$mean_eCDF_diff[dat$X==3], dat$mean_eCDF_diff[dat$X==16]), type="l", lwd=.8, col="darkgray", pch=20, cex=.5) 
lines(c(1.03,1.97), c(dat$mean_eCDF_diff[dat$X==4], dat$mean_eCDF_diff[dat$X==17]), type="l", lwd=.8, col="darkgray", pch=20, cex=.5) 
lines(c(1.03,1.97), c(dat$mean_eCDF_diff[dat$X==5], dat$mean_eCDF_diff[dat$X==18]), type="l", lwd=.8, col="darkgray", pch=20, cex=.5) 
lines(c(1.03,1.97), c(dat$mean_eCDF_diff[dat$X==6], dat$mean_eCDF_diff[dat$X==19]), type="l", lwd=.8, col="darkgray", pch=20, cex=.5) 
lines(c(1.03,1.97), c(dat$mean_eCDF_diff[dat$X==7], dat$mean_eCDF_diff[dat$X==20]), type="l", lwd=.8, col="darkgray", pch=20, cex=.5) 
lines(c(1.03,1.97), c(dat$mean_eCDF_diff[dat$X==8], dat$mean_eCDF_diff[dat$X==21]), type="l", lwd=.8, col="darkgray", pch=20, cex=.5) 
lines(c(1.03,1.97), c(dat$mean_eCDF_diff[dat$X==9], dat$mean_eCDF_diff[dat$X==22]), type="l", lwd=.8, col="darkgray", pch=20, cex=.5) 
lines(c(1.03,1.97), c(dat$mean_eCDF_diff[dat$X==10], dat$mean_eCDF_diff[dat$X==23]), type="l", lwd=.8, col="darkgray", pch=20, cex=.5) 
lines(c(1.03,1.97), c(dat$mean_eCDF_diff[dat$X==11], dat$mean_eCDF_diff[dat$X==24]), type="l", lwd=.8, col="darkgray", pch=20, cex=.5) 
lines(c(1.03,1.97), c(dat$mean_eCDF_diff[dat$X==12], dat$mean_eCDF_diff[dat$X==25]), type="l", lwd=.8, col="darkgray", pch=20, cex=.5) 
lines(c(1.03,1.97), c(dat$mean_eCDF_diff[dat$X==13], dat$mean_eCDF_diff[dat$X==26]), type="l", lwd=.8, col="darkgray", pch=20, cex=.5)  
text(0.74,dat$mean_eCDF_diff[dat$X==1],"propensity", font=1, cex=0.8) 
text(0.69,dat$mean_eCDF_diff[dat$X==2]+0.01,"serious fraud", font=1, cex=0.8) 
text(0.64,dat$mean_eCDF_diff[dat$X==3]-0.005,"loser's vote share", font=1, cex=0.8) 
text(0.79,dat$mean_eCDF_diff[dat$X==4]-0.002,"boycott", font=1, cex=0.8) 
text(0.55,dat$mean_eCDF_diff[dat$X==5]-0.01,"poll type, leader tenure", font=1, cex=0.8)  
text(0.56,dat$mean_eCDF_diff[dat$X==6]-0.01,"pre-election violence", font=1, cex=0.8) 
text(0.70,dat$mean_eCDF_diff[dat$X==8]-0.02,"post-conflict", font=1, cex=0.8) 
text(0.44,dat$mean_eCDF_diff[dat$X==9]+0.01,"GDP pc, executive constraints", font=1, cex=0.8) 
text(0.60,dat$mean_eCDF_diff[dat$X==10]-0.01,"natural resources pc", font=1, cex=0.8) 
text(0.67,dat$mean_eCDF_diff[dat$X==12]+0.01,"population size", font=1, cex=0.8) 
text(0.54,dat$mean_eCDF_diff[dat$X==13]+0.01,"ethnic fractionalization", font=1, cex=0.8) 
legend(0.95, 0.8, bty="n", c("difference significant","difference not significant"), 
       col=c("black","darkgray"), cex=0.9, pch=c(19,19))
dev.off()



#######################
# Figure A4
#######################
# code adapted from Matt Pietryka, which is adapted from VanderWeele 2011
rm(list=ls(all=TRUE))
library(foreign)
library(readstata13)
library(rms)
dat<-read.dta13("Condemnation.dta")
attach(dat)

mod2=lrm(PostElViolence ~ NegReport + observedW46 + SeriousFraud + LoserVoteShare + Boycott + type3 + previolev + DPI_yrsoffc_Lag1Log + Prio_CWPast10 + WB_GDPpc_Lag1Log + Ross_oilgas_valuePOP2000_Lag1Log + P_xconst_Lag1 + WB_PopM_Lag1Log + Finkel_EF_Lag1, x=T, y=T, data = dat)
mod2
robcov(mod2, cluster=dat$cowcode)

CondemnCoef <-2.8786
CondemnSE <- 0.8925
CritVal <- 1.96   # approx 95% CI
CondemnU <- CondemnCoef + (1.96 * CondemnSE)
CondemnL <- CondemnCoef - (1.96 * CondemnSE)
Estimates<-cbind(CondemnU, CondemnCoef, CondemnL)
print(Estimates)

# GLM results
s.glm <- function(the.coef,                                 #coef of interest 
                  gamma.pos=c(1,1.5,2,2.5,3,3.5,4, 4.5,5),  #effect of latent factor on Y
                  pi.1=.9,                                  #prevalence of latent factor when X = 1
                  pi.0=.1                                   #prevalence of latent factor when X = 0
) {
  if(the.coef>=0){
    ###SENSITIVITY
    odds.ratio<- exp(the.coef)                              #the estimated effect, expressed as odds ratio
    the.bias<-(1+(gamma.pos-1)*pi.1)/(1+(gamma.pos-1)*pi.0) #the bias (when beta is positive)
    corrected.n.t<-odds.ratio/the.bias ;round(corrected.n.t,2) #the estimate corrected for assumed level of bias
    return(rbind(gamma.pos,corrected.n.t))
  }
  else{
    gamma.neg<-1/gamma.pos                                  #the odds ratio of the effect of unmeasured factor on Y (when beta is negative)
    ###SENSITIVITY
    odds.ratio<- exp(the.coef)                              #the estimated effect, expressed as odds ratio
    the.bias<-(1+(gamma.neg-1)*pi.1)/(1+(gamma.neg-1)*pi.0) #the bias (when beta is negative)
    corrected.n.t<-odds.ratio/the.bias ;round(corrected.n.t,2) #the estimate corrected for assumed level of bias
    return(data.frame(gamma.neg,corrected.n.t))
  }
}

s.glm (CondemnCoef)
s.glm (CondemnU)
s.glm (CondemnL)

Corrected<-rbind(s.glm (CondemnCoef), s.glm (CondemnU), s.glm (CondemnL))
transposed<-apply(Corrected,1,t) 
colnames(transposed) <- c("Gamma","CorrectedCoef","","CorrectedU","","CorrectedL")
Estimates<-as.data.frame(transposed)

# Replicate analysis for only OBSERVED elections 
observedEl <- dat[ which(dat$observedW46 == 1 ), ]
mod3=lrm(PostElViolence ~ NegReport + SeriousFraud + LoserVoteShare + Boycott + type3 + previolev + DPI_yrsoffc_Lag1Log + Prio_CWPast10 + WB_GDPpc_Lag1Log + Ross_oilgas_valuePOP2000_Lag1Log + P_xconst_Lag1 + WB_PopM_Lag1Log + Finkel_EF_Lag1, x=T, y=T, data = observedEl)
mod3
robcov(mod3, cluster=observedEl$cowcode)

CondemnCoefObs <-3.5454 
CondemnSEObs <- 1.2964
CritVal <- 1.96   # approx 95% CI
CondemnUObs <- CondemnCoefObs + (1.96 * CondemnSE)
CondemnLObs <- CondemnCoefObs - (1.96 * CondemnSE)
EstimatesObs<-cbind(CondemnUObs, CondemnCoefObs, CondemnLObs)
print(EstimatesObs)

s.glm (CondemnCoefObs)
s.glm (CondemnUObs)
s.glm (CondemnLObs)

CorrectedObs<-rbind(s.glm (CondemnCoefObs), s.glm (CondemnUObs), s.glm (CondemnLObs))
transposedObs<-apply(CorrectedObs,1,t) 
colnames(transposedObs) <- c("Gamma","CorrectedCoef","","CorrectedU","","CorrectedL")
EstimatesObs<-as.data.frame(transposedObs)

# Figure A4 
pdf(file = "FigureA4.pdf",
    width = 8, height = 8, onefile = TRUE, family = "Times")
par(mar=c(4,4,2,2)+.1)

plot(EstimatesObs$Gamma, EstimatesObs$CorrectedCoef, type="b", lty=1, pch=23, col="black", ylim=c(0,35), xlab="Odds Ratio for Omitted Variable (Gamma)", ylab="Corrected Odds Ratio for Condemnation", yaxt="n")
axis(2, las=1)
lines(EstimatesObs$Gamma, EstimatesObs$CorrectedL, type="b", lty=3, "black", pch=23)
lines(Estimates$Gamma, Estimates$CorrectedL, type="b", lty=3, pch=17, col="blue")
lines(Estimates$Gamma, Estimates$CorrectedCoef, type="b", lty=1, pch=17, col="blue")
legend(2.1, 35, bty="n", c("Estimate (Observed Elections Only)","Lower Bound (Observed Elections Only)",
                           "Estimate (All Elections)","Lower Bound (All Elections)"), 
       col=c("black","black","blue","blue"), cex=1.2, lty=c(1,3,1,3), pch=c(23,23,17,17))
dev.off()
